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Many real oscillators are coupled to other oscillators and the coupling can affect the response of 
the oscillators to stimuli. We investigate phase response curves (PRCs) of coupled oscillators. The 
PRCs for two weakly coupled phase-locked oscillators are analytically obtained in terms of the PRC 
for uncoupled oscillators and the coupling function of the system. Through simulation and analytic 
methods, the PRCs for globally coupled oscillators are also discussed. 
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Many systems in physics, chemistry and biology are 
modeled as interacting nonlinear oscillators [J, 0, y|, 0, H, 
6] . One of the easiest ways to characterize an oscillator is 
its phase response curve (PRC) 0, 0, H, H, 0] ■ The PRC is 
defined as the steady phase shift of an oscillation relative 
to the unperturbed oscillation as a function of the timing 
of perturbation to the oscillator. It provides a useful 
information for understanding the oscillator's behavior 
when the oscillator is subjected to external stimuli or 
signals from other oscillators. 

In most of previous studies, the PRC is obtained when 
the oscillator is isolated from other oscillators [H, 0, 0, 0] • 
However, many oscillators in real systems are coupled to 
others when they are under the influence of external stim- 
uli, and the coupling can affect the response of the oscil- 
lators. To better understand the dynamics of oscillators 
such as the response of neuronal population to signals 
from other brain region Q or to controlling stimulations 
Q, it is necessary to study how the coupling changes 
the PRCs. This study can also give insights into the 
phase response of a giant oscillator (for example, circa- 
dian rhythm generators Q ) composed of many individual 
oscillators Q. In this letter, we study the PRC of cou- 
pled oscillators using the average phase of the system and 
the relative phases between the oscillators comprising the 
system. The PRC is shown to depend on the PRC of the 
isolated oscillator, the nature of the coupling, and the 
relative phases between the oscillators. For some cases, 
the PRCs are analytically obtained. Our approach differs 
from that of Ref. [8| in that we analytically approximate 
the PRC while they require the numerical evaluation of 
the adjoint of a certain linear operator. 

If coupling between a network of oscillators is suffi- 
ciently "weak" , the possibly high-dimensional system can 
be reduced to a network of coupled phase models 0, 0, 0] • 
In the following wc exploit this fact and restrict our anal- 
ysis to coupled phase models. Consider, first, two weakly 
coupled phase-locked oscillators subjected to a common 



perturbation characterized by their individual PRC: 

0i = w 1 + KH{e 2 -6 1 ) + Z(9 1 )-A6(t-t 1 ), (1) 
6 2 = uj 2 + KH(9 1 - 6 2 ) + Z{6 2 ) ■ AS(t - t x ), (2) 

where 6i(t) is the phase of oscillator i at time t, u>i is the 
natural frequency of the oscillator i and K(> 0) is the 
coupling strength. H{9) is the coupling function obtained 
by the phase reduction [1, 0, 0]. A5(t — ti) denotes a 
Dirac delta impulse with amplitude A at time t\ which is 
sufficiently large so that the perturbing impulse is applied 
after the system reaches a steady state. Z{&) is the PRC 
for uncoupled oscillator obtained using an impulse with 
unit amplitude. Without coupling (K = 0), the impulse 
causes steady phase shift AZ{0i(t\)) for oscillator i. 

In the presence of coupling (K ^ 0), if the oscillators 
are locked with nonzero phase difference, or the input 
amplitudes are different, then the input impulse gener- 
ally causes nonidentical phase changes to the oscillators. 
Thus, the system transiently deviates from the locked 
state and then returns to the state. The coupling can af- 
fect the phase shift which the oscillation of the recovered 
state can have relative to the unperturbed oscillation. Wc 
wish to determine the PRC of the coupled oscillators, in 
other words, how the phase shift depends on the phase 
at ti of the perturbation. 

To analyze the dynamics, we convert Eqs. (|TJ) and 
@ into those for the average phase $ s 
relative phase (f> = 6\ — 2 . 
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$ = uo + KH e (^) + Z av {0 u 6 2 )- ASit-tr), (3) 
a\ = Auj-2KH o (<p) + Z d (0 1 ,8 2 )-A5(t~t 1 ), (4) 



where u> = 
z(e 1 )+z(e 2 ) 
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For simplicity, let us assume that the system has one 
stable locked state with cf) = <f>Q satisfying = Aoj — 
2KH (<j>o) and H '{4>q) > 0. The phase of each oscillator 
can be written as 6\ = $ + (f>/2 and 2 = $ — <fi/2. Let us 
denote the phase shift in a phase, for example 6\, relative 
to the unperturbed oscillation by A6*i. We can see that 
the phase shift for the oscillator 1 is given by 
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A0i = A$ + A0/2. 



(5) 
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FIG. 1: (Color online) PRCs with odd coupling functions. 
lui = tv/2 + Auj and co 2 = tt/2. Auj = 0.4 and A = 1. 
(a) K = 1: 00 ~ 0.201. (b) K = 0.25: O « 0.927. (c) 
A' = 0.22: 0o ~ 1.141. For (a), (b), and (c), H{6) = sinfl 
and Z(6>) = -sin<9. (d) #((9) = sin (9 - 0.4sin(2<9), Z{6) = 
- [sin(0 + 0.27T) -sin(0.27r)]/[l + sin(0.27r)], and K = 0.5: 
0o ~ 0.908. 

When H{9) is an odd function, the average phase <& 
evolves with a constant frequency uj before and after the 
impulse (Eq. ©). Thus, A$' = AZ av {6 l {t 1 ),6 2 (t l )). 
When the relative phase remains in the basin of attrac- 
tion of the original relative phase right after the impulse, 
approaches the original relative phase. Otherwise, the 
relative phase moves to another stable value (called walk- 
through). Thus, A0 = <fif — 0o, where 0/ is the stable 
value of reached after the impulse. Note that even 
= 0o and = 0o + 2ir give different results. Therefore, 
the PRC of the oscillator 1 in the coupled cases is given 
by 

Zd(0i) = AZ av (9 1 ,e 1 - O ) + (0/ - o )/2. (6) 

We simulate Eqs. ([TJ) and @ using Euler method with 
time step At = 0.01. We measure the steady phase shift 
due to the impulse relative to the unperturbed activity. 
The PRC is given by this phase shift as a function of the 
phase at which the impulse is applied. 

Figure Q] shows Z C \{Q\) with odd coupling functions. 
The prediction from the theory (black solid curves) 
matches very well with the simulation results (symbols). 
With larger values of Aw and/or smaller values of K, the 
oscillators are locked with larger 0o. In Figs. [TJa), (b), 
and (c) with H(9) = sin# and Z(9) = — sin(9, we show 
the PRC for different values of the coupling strength K. 
When 0o is very small, the PRC of coupled oscillators 
is very close to that of uncoupled oscillators as expected 
(Fig. Ha)). In this case, goes to the original value 
0o after the impulse. In Fig. HJb), with the larger 0o, 
the PRC of the coupled oscillator becomes significantly 
different from that of uncoupled oscillators. When the 
impulse can kick the system out of the basin of the sta- 
ble locked state with 0o, the system goes through phase 
walk through. If the system has a stable fixed point with 



FIG. 2: (Color online) PRCs with non-odd coupling 
functions. For (a)-(c), u\ — n/2 + Auj, u>2 = n/2 
and Auj = 0.4. (a) H{6) = sin(6> + 0.4tt), Z{6) = 
-[sin(6» + 0.47r) -sin(0.47r)]/[l + sin(0.47r)], and K = 1: 
0o ~ 0.704, A = 1. (b) H{6) = sin(fl - 0.4?r) + 0.3sin(26> - 
O.Itt), Z(6) = -[sm(0 + O.27r)-sin(O.27r)]/[l + sin(O.27r)], 
and K = 0.4: O ~ 0.769, A = 0.2. (c) H{6), Z(0), and 
K are the same as in (b) : 0o ~ 0.769, A — 0.5. (d) gap- 
junction coupled Morris-Lecar oscillators Qj: I^ x t,i = Io + AI 
and hxt,2 = h - AI with AI = 0.2, A v = 40. (top) type I 
case, Io = 50 and g 3yn = 0.015, 0o ~ 1.279. (bottom) type II 
case, Io = 94 and g syn = 0.01, 0o « 0.778. 

00 and an unstable fixed point 0„ in [0, 27r) as in the 
case with H(6) = sin# for Aw < 2K, <fi u has the role 
of basin boundary and goes to 0/ = 0o + 27r when 
00 + AZ (0i(h)) -AZ(6 2 (ti)) >0„. This type of changes 
in causes the discontinuity shown in the PRC of Fig. 
[IJc). We show similar results for a coupling function with 
higher order Fourier terms and an asymmetric PRC (Fig. 

When H(9) is not an odd function, the even part of 
H affects the dynamics of $ and thus the phase shift 
A$ through Eq. ©. Finding the PRC in the ana- 
lytic form is not possible for these cases since we have 
to solve equation for general initial data. Instead, 
we can get an approximation of the PRC in the limit 
of small changes in 0. Let = 0o + q with \q\ <C 1. 
We can linearize Eq. ([4]) and obtain approximation 
q « qoe~ 2KH ° (^°K i-tl ) for t > t\ where qa is the change 
in right after the impulse: q n = AZd(9i(ti),82(ti)). 
As t — > oo, returns to O . Thus, A0 = 0. The 
phase shift A5> is given by A$ = AZ av (6 1 (t 1 ),8 2 (t 1 )) + 
f™K[H e (ct>)-H e ((t> )}dt w AZ av (9 1 (t 1 ),e 2 (t 1 )) + 

2H o '(M q0 > Where Wt3 USe He ^ ~ ff eOo) ~ H e '((j) Q )q(t). 

Therefore, the PRC of the oscillator 1 is 
Z cl {6 x ) w AZ av {9x,Bi-<h) 

2H (0 O ) 

Figure [2] shows Z c i{8{) with non-odd coupling func- 
tions. In Fig. [DJa), we show the PRC with the sim- 
ple type of H. While AZ and AZ av are similar, the 
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obtained PRC for the coupled oscillator is significantly 
different from them. The curve from Eq. ([7]) fits well 
with simulation results for the entire range of 9\. Fig- 
ure [DJb) shows the results with a H function with higher 
order terms. We use small A = 0.2 for this case and 
the PRC from the theory fits well with the simulation 
result. In Fig. EIc), we use the same parameters as 
in (b) except A = 0.5. With the larger A, the theory 
mismatches significantly for a range of phases, but still 
gives a relatively similar shape to the simulations. The 
overall matching is due to the fact that q = at some 
phases satisfying Z(9) = Z(6 — </>o) and around those 
phases the theory fits well with simulation results. Fig- 
ure Efd) shows the PRC of gap-junction coupled Morris- 
Lecar oscillators with slightly different injection currents 
i: CVi = -IiV^wJ+I^i+g^niVj-ty+AySit-h) 
with j — 2,1 for i = 1,2. The details are in Ref. 

The system is simulated using the 4th-order Rungc- 
Kutta method. For type I and type II cases, the theory 
gives good fitting with the simulation results with weak 
stimulus. 

Next, we want to understand PRCs for oscillators cou- 
pled to many other oscillators. We study the case with 
globally coupled oscillators: For i = 1,2, N, 

K N 

k = "i + jjYl H & - *<) + Z ^ ■ M ^ - (8) 

i=i 

where N is the total number of oscillators and others are 
as defined in the two oscillator system. 

We introduce similar variables as in two coupled oscil- 
lators: the average of the phases $ = ^2j=i ®j an0 - ^ ne 
phase 4>i = 6i — 9m of oscillator i relative to the phase 
of oscillator M, where the subscript M denotes the os- 
cillator which has the average frequency Co. From the 
definitions of $ and 4>, we obtain 9m = < E > — Ylj=i §j- 
Because the PRC for other oscillators can be treated sim- 
ilarly and oscillator M follows closely to the collective 
behavior of the system, we focus on the PRC Z c (9m) of 
oscillator M. As in the case of two coupled oscillators, 
we get 

1 N 

Z C {9 M ) = A$ - (A<j>) with (A<t>) = jrJ2 A fa- W 

3=1 

The equation for <E> is 
K N 

$ = Q + — ^2 H e (<t>j - 0j) + Z av ■ AS(t - h), (10) 

i,j=i 

where Z av {9 x , ...,9 N ) = ± ^ti Z{9 l (t))- 

For an odd function H(9), A$ = AZ av . But for a 
non-odd function H(9), the second term contributes to 
A$ and it is not easy to calculate A$ analytically. 

Let us consider fully locked states first. For a fully 
locked state with relative phases fyo, the system returns 
to the locked state after the stimulation and the relative 



phase 4>iQ can be changed to the equivalent phase <piQ + 
2riiTr, where rii is an integer. Thus, Afa = 2nj7r. 

With H{9) = sin(0+/3), which is a good approximation 
for many general coupling functions, Eq. ([5]) becomes 

6\ = co t + KRsm(e - 9 t +/3) + Z(9 t ) ■ A5(t - t x ), (11) 

where R and are the order parameter and the corre- 
sponding collective phase respectively defined by Re = 
■h Y^jLi e%8i ■ I n the frame rotating with the synchroniza- 
tion frequency SI, the equation becomes 

fa = oj 1 - SI + KRs'm(Q - fa + (3) 

+Z(6 i )-A5(t-t 1 ), (12) 

where fa = 6i - (Sit + O ), = 6 - (Sit + O ) and the 
constant Oo is chosen such that the stationary value of 
before the impulse is zero. Let Rq denote the stationary 
value of R. We can analyze the stationary state of the 
system, using self-consistency argument and find Rq and 
Si 

For fully locked states or partially locked states where 
the oscillator M is locked with a locking phase faw* and 
the oscillators form a stationary distribution relative to 
the frame rotating with SI, Z av = Xa=i Z(9m + ipi — 
ip M *), where Va/* = sin^f^) + (3. For Z(0) = oi + 
a 2 sm(6* + £), using R = X)j=i e *' we obtain 

Zav = ai + a 2 R sm(9M + £ - i'M*)- (13) 

Note that since Rq < 1, the magnitude of the sinusoidal 
part of Z av is smaller than that of Z unless synchrony is 
perfect. 

Figures Oa)-(c) show results with H(9) = sin 9 and a 
uniform distribution for the frequencies of the oscillators. 
We use Z(6) = — sin (9 for (a)-(d), and an asymmetric 
Z(6) for (e). With the given coupling strength, the sys- 
tem shows a fully locked state (Fig. [3ja)). Figures EJb) 
and (c) show the PRCs for different values of A. With 
weak stimulation (Fig. [3jb)), all <f>i return to the unper- 
turbed values (Afa = for all i, Fig. [21(d)) and the PRC 
is shown to be contributed only by A<f> = AZ av . The 
prediction Z c (0m) = AZ av from the theory (Eq. ([13])) 
fits well with the simulation results. In contrast, with a 
stronger impulse (Fig. Efc)), the simulation results de- 
viate from Z c (9m) = AZ av for some range of 9m- We 
calculate (A<p) from the simulations and it accounts for 
the deviation as predicted from Eq. ©. 

The deviation in Fig. [3](c) can be understood as fol- 
lows. Nonzero Afa can occur only when the order pa- 
rameter transiently decreases. For 9m G (7r/2,vr), the 
impulse disperses the locked group, because the trailing 
oscillators receive more negative impact than the lead- 
ing ones. Thus, the order parameter decreases from Rq 
to R(ti+) (< Rq). R can decrease more depending on 
the behavior of the oscillators, and then returns to i?o- 
The collective phase O also decreases (0(ti+) < 0), be- 
cause most of the phases of the oscillators decrease due 
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FIG. 3: (Color online) Fully locked cases, (a) snapshot of 
phases of oscillators: a fully locked state (b) A — 0.5. (c) 
A = 1.0. (d) Afa for (b) and (c) with 6 M = 0.825tt. For 
(a)-(d), H(9) = sine and Z(9) = -sin0. (e) H(9) = 
sin 9, Z(9) = - [sin(0 + 0.3vr) - sin(0.3ir)] / [1 + sin(0.37r)] , 
and A = 0.5. (f) H{9) = sin(6» + O.Itt), Z(9) = -sin6>, 
and A = 0.5. A uniform distribution is used for {c^i}: (i) 
uii — lj — Au + 2A ^ < ^~ 1 ^ or (ii) randomly selected ua from 
[uj — Alj, uj + Au>] . (i) is used for the symbols of (a)-(e) except 
the gray circles in (b), (c), (e), and (f). U) = tt/2. Auj = 0.6. 
K = 0.8 for (a)-(e) and K = 1 for (f). 

to the impulse. The sudden changes in R and af- 
fect the dynamics of oscillators. The behaviors of oscil- 
lators right after the impulse can be described by the 
equation ^ = — il + KR(tx+) sin(0(ti+) — ipi) with 
ipi(h+) = ipi* + Z(6 M (h) ~i>M*), where ipi* is the 
locking phase for oscillator i. The trajectory of oscillator 
i can escape completely from the basin of attraction of 
fao during the transient behavior of R and settle to the 
equivalent phases fa = fao + 2niir. Since the curves for 
(ipi, tpi) arc shifted to the left due to the negative 0(f i+) 
and upwards (downwards) for the oscillators with u>i > 
(u>i < 17), the oscillators with frequencies far from the 
average one can escape and those with higher frequen- 
cies escape first. Because of this, (Afa > and the PRC 
deviates negatively from AZ av (Eq. (flU)) ). The oscilla- 
tors with frequencies far from the average one have more 
chance to have higher n (Fig. [3]Jd)), because they can 
drift faster and stay unlocked longer. Other ranges of 
8m can be understood similarly. 

When H{9) is not odd, it is difficult to find the gen- 
eral results. For H(0) = sin(# + (3), we can see that the 
second term of Eq. (flT))) is equal to KR 2 sinf3. Thus, 
A$ = AZ av + K sin j3 f™(R 2 - R 2 )dt. With weak 



stimulus, Z c (6m) = A$ and the PRC deviates posi- 
tively (negatively) from AZ av for values of R(ti+) > Rq 
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FIG. 4: (Color online) Partially locked cases with H(9) = 
sin 9. (a) snapshot of phases of oscillators: a partially locked 
state, (b) A = 0.5. (c) Afa. Z(9) = -sin0. {uJi} 

obeys a Gaussian distribution g(u), Q) = J 2 exp(— ^" 2 ~"^ ): 
(i) W(jv/2)+* = Q + yk, W(jv/2)-fc+i = & - Vk with y k = 
(xh-i + x k )/2, x k +i = Xk + N~ 1 /g(x k , 0), and x = for 
k = 1, N/2 or (ii) = Q - yt and Lu i+N / 2 = u> + yi 
with i < N/2 and yi randomly selected according to g(y,0) 
with y > 0. (i) is used for the simulations(symbols) of (a)-(c) 
except the gray circles in (b). Gj = n/2, a — 0.3, and K = 0.6. 



(R(h+) < R ) (Fig. Etf)). The values of R(h+) arc 
easily calculable using Z and the distribution for uj. 

Finally, let us briefly consider partially locked cases in 
the limit of — ► oo. When the system exhibits a par- 
tially locked state, the drifting oscillators form a station- 
ary distribution in the frame rotating with the synchro- 
nization Q. In the original frame, we can say that the dis- 
tribution rotates with f2. We can define PRCs for locked 
oscillators. Let us consider cases with H(9) = sin (9. Fig- 
ure 0Jb) shows the PRC Z C (9 M ) of oscillator M for the 
partially locked state of (a). We can understand Z c (0m) 
through Eq. ([9]). Since H is an odd function, A$ = AZ av 
(Eq. ([TH]) ). While for the locked oscillators Afa = In^ 
and is nonzero in some ranges as in the fully locked cases, 
for the drifting oscillators is usually not an integer 
multiple of 2tt and can be nonzero in any ranges (Fig. 
0Jc)). Simulations show that Afa for the drifting oscil- 
lators contribute significantly to the PRC and the PRC 
(symbols, Fig. 0Jb)) differ from AZ av (the dashed curve) 
for almost the entire range of 8m- 

In summary, we have investigated the PRCs of coupled 
oscillators in terms of the PRCs of individuals, the nature 
of the coupling, and the relative phases of the oscillators. 
Our approach of obtaining PRCs using the average and 
relative phases can be applicable to oscillators on differ- 
ent type of networks. 

This work was supported by National Science Founda- 
tion grant DMS05135. 
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